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NUMERICAL STUDY OF SOUND PROPAGATION IN A JET FLOW 


Sharon L. Padula and Chen-Huei Liu* 

Langley Research Center 

SUMMARY 

An efficient and reliable numerical method for solving problems related to sound 
propagation through a jet flow field is demonstrated by applying the method to a simple 
test problem. The test problem is reformulated in terms of a slowly changing variable 
and a finite -difference approach is employed to solve the resulting elliptic-type partial 
differential equation. This discretization process yields a system of linear algebraic 
equations which can be solved by a standard library subroutine with decreased computer 
time and storage requirements. Comparison of numerical results with those of previous 
studies and with experimental data demonstrates the improved accuracy of the present 
analysis. 


INTRODUCTION 

Understanding the principles of noise propagation is an essential ingredient of sys- 
tematic jet noise reduction research. High-speed computer methods offer a unique poten- 
tial for dealing with complex real life physical systems whereas analytical solutions are 
restricted to simplified idealized models. The formulation of sound propagation through 
a jet flow in terms of pressure was found to be impractical for computer solutions and a 
more suitable approach was needed. Previous investigations by Schubert (ref. 1) and by 
Liu and Maestrello (ref. 2) selected the phase and amplitude of the acoustic pressure as 
dependent variables requiring the solution of a system of nonlinear algebraic equations. 

The nonlinearities complicated both the analysis and the computation. A reformulation of 
the convective wave equation in terms of a new set of dependent variables was suggested 
(ref. 3) with a special emphasis on its suitability for numerical solutions on fast computers. 
The technique is very attractive because the resulting equations are linear in a slowly 
changing variable. The computer solution to such a linear system of algebraic equations 
may be obtained by well defined and direct means which conserve computer time and stor- 
age space. 

In this paper, the reformulation technique is applied to cases equivalent to those pre- 
sented in reference 1. Complete details of the formulation are contained in the appendix. 

*NRC-NASA Resident Research Associate. 



Numerical results from the present analysis and from reference 1 are compared with 
applicable experimental measurements by Grande tref. 4) and with similar numerical 
results of Mungur et al. (ref. 5). 

SYMBOLS 

A amplitude of radiated sound wave 

c ambient speed of sound 

C nondimensional c, c/uj 

d jet exit diameter, 0.01905 m 

f frequency 

H, K step sizes, and At], respectively 

M local flow Mach number, u/c 

Mj flow Mach number at jet exit, u-j/c 

n,m number of grid lines 

p pressure variable 

r ,9,<p spherical polar coordinates (see fig. 1) 

R nondimensional r, r/d 

t time 

T nondimensional t, tuj/d 

u longitudinal component of flow velocity 

uj jet exit velocity 

U nondimensional u, u/uj 

x Cartesian coordinate in windward direction of jet 
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X nondimensional x, x/d 

Z complex variable, AR exp(i^j) 

77 coordinate transformation, tan - * ( a ~ 2 tan Q) 

£, coordinate transformation, In R 

a parameter of angular stretching, 0.538 (see eqs. (5)) 

$ velocity- potential variable 

xp phase of radiated sound waves, \pQ + xpi/R 

components of xp 
co angular frequency 

£2 nondimensional angular frequency, cud/uj 

FORMULATION OF PROBLEM 

For illustration purposes, the method developed in this paper is applied to a simple 
test problem. In a subsequent section, the appropriate convective wave equation is identi- 
fied and is reformulated for solution. A judicious choice of variables nets a linear partial 
differential equation which can be efficiently solved by numerical methods. 

The Nonuniform Jet Flow Field 

The present method of solution is tested for the case of an acoustic point source 
located on the center line and contained within the potential core of a spreading jet. 

(See fig. 1.) The jet nozzle is cylindrical with diameter d. The jet flow is assumed 



Figure 1.- Jet flow configuration. 
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unidirectional and axisymmetric but nonuniform and the velocity profiles at any cross 
section are determined from experimental and empirical data. The origin of the coordi- 
nate system is fixed at the source while the origin of the flow-field coordinates is fixed 
at the virtual origin of the jet. Most of the calculations are done in the spherical polar 
coordinates (r ,0, <£). 


The Convective Wave Equation 

Sound propagation through a flow field can be described most simply by the following 
convective wave equation: 


L(*) = 0 


( 1 ) 


where 


L = 2u-g- 2 - 

c 2 \^2 dx dt 


+ u* 


- V* 


is a linear differential operator, $ is Obukhov’s "quasi-potential" variable, v2 is the 
Laplacian operator in Cartesian coordinates, c is the ambient speed of sound, and u is 
the longitudinal component of flow velocity. A discussion of the use and validity of this 
formulation can be found in reference 1. 

By assuming a harmonic source, $ can be expressed as 


$ = A exp [i(i//R - OT)] 


( 2 ) 


where £2 = cod/uj is the normalized source frequency, R*r/d is the radial distance 
from the source, T = tuj/d is the time, and A and if/ are the amplitude and phase of 
the radiated sound waves. 

Substituting $ into equation (1) results in a time -independent partial differential 
equation which is nonlinear in A and if/. Schubert (ref. 1) chooses to solve this equation 
by application of the finite -difference method. This necessitates solving a system of simul- 
taneous nonlinear algebraic equations. 

There is no direct method for solving such a system of nonlinear equations. Iterative 
computer methods are available, but they have a number of drawbacks. First, they are 
highly specific. The coding is dependent on the problem to be solved and must be rewrit- 
ten and tested for each new equation. Secondly, the accuracy of the results can be affected 
by how well and how quickly the process converges. In the case of sound propagation 
through a flow, it remains to be shown that the process does in fact converge for all 
instances. 
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The problems inherent in iterative methods can be avoided by representing the 
variable i p by 

^ = ^0 + ( 3 ) 

where is a function of r and 9 and \j/ approaches i//q as R approaches infin- 

ity. Then substituting equation (3) into equation (2) gives 

$ = |exp[i(^°R - f2T)j (4) 

where Z is a complex variable equal to AR exp(ii^. Equation (1) can be rewritten as 
a linear partial differential equation in Z. (See eq. (A2) in the appendix.) In this way, the 
solution by finite -difference method involves merely the solution of a system of simulta- 
neous linear algebraic equations. Most computer subroutine libraries contain routines 
adequately suited to this task. 


NUMERICAL SOLUTION 

In this section, a solution to the test problem by numerical means is discussed. The 
techniques employed help to reduce further the computer storage space and central mem- 
ory time required. 


Underlying Assumptions 

The basic assumptions and basic scheme for solving the problem follow directly from 
previous work. (See refs. 1 and 2.) Details of the formulation can be found in the appen- 
dix. Figure 2 shows the spreading jet superimposed on a polar grid. An antijet is assumed 
in order to avoid problematic boundary conditions along the rigid walls of the jet nozzle. 

The value of the complex variable Z changes most rapidly near the point source 
and in the region of significant flow. It is desirable to apply the finite-difference method 
to an unequally spaced grid which has a concentration of points in the areas of greatest 
change. (See fig. 2.) To avoid the difficulties posed by an uneven grid, an even grid is 
specified in some new coordinates, £ and t). (See fig. 3.) This new grid in (£, 77 ) maps 
onto the desired uneven grid in (R, 9) according to the following transformation: 


| = In R 

7] = tan"* (a - 2 tan 9) 


(5) 
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Here, a is a parameter (<j 2 < i) which determines the amount of angular stretching; 
a = 0.538 throughout this paper. For illustration purposes, the polar grid in (£,rj) is 
thought of as a rectangular grid with n x m intersecting lines. 

The Finite-Difference Method 

The partial differential equation (1) written in terms of Z should be satisfied at 
each interior point of the region of interest. If the partial derivatives in this equation are 
replaced by their central difference approximations, a nine -point difference equation 
results. Writing this difference equation at each interior grid point gives (n - 2) x (m - 2) 
equations in n x m unknowns. Boundary conditions are the additional equations needed 
to specify a unique solution. 


Boundary Conditions 

The equation in Z is an elliptic-type partial differential equation, thus, the 
boundary -value problem is well posed. 

The region T) is bounded by two concentric half circles centered on the point 
source and by two line segments on the axis of symmetry of the jet. (See fig. 4.) The 
inner circle has a radius of 0.25d and the outer circle has a radius of lOOd, where d is 
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the diameter of the jet nozzle. In this way, the inner boundary is contained in the potential 
core of the jet, and the outer boundary approximates the far field. 

The inner boundary conditions come from a solution derived by Moretti and Slutsky 
(ref. 6). The solution is specified for a point source in a uniform flow. It gives the sound 
pressure level for points a small distance from the source. Since the acoustic source and 
the inner boundary are contained in the potential core of the jet, they are in a locally uni- 
form flow and the Moretti and Slutsky solution is approximately valid. 

The outer boundary conditions at lOOd are essentially the far-field radiation condi- 
tions. Along the center line, symmetry conditions are imposed in discrete form. 
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RESULTS AND DISCUSSION 


Numerical results produced by the present analysis are presented in figures 5 to 11. 
Test cases which demonstrate the utility and reliability of the analysis for a wide range of 
flow Mach numbers, source frequencies, and radial distances are presented first. The 
graphs contain either sound pressure level (SPL) or phase information. Comparisons with 
available experimental data and numerical results follow. 

Test Cases of the Present Analysis 

Unless otherwise stated, all test results have been calculated on a grid of 
19 x 19 lines, i.e., n = m = 19. Figure 5 indicates that for the problem tested a denser 
grid spacing will cause very little change in the results and reasonable approximate results 
can be achieved with an 11 x 11 or 13 x 13 grid. Such results can be produced very rapidly 
in a minimal amount of core storage space. 



Figure 5.- Effects of the number of grid lines on the solution. 


Figure 6 is representative of the numerical results possible with the current 
approach. Constant sound-pressure-level curves are plotted for a case in which the 
acoustic point source is located at the intersection of the jet exit plane and the jet center 
line. Here the flow velocity was Mach 0.3 and the frequency was 3000 Hz. 

For the purpose of comparing phase information produced by the present method 
with that produced by reference 1, results over a range of flow Mach numbers and at three 
different radial distances from the source are presented in figure 7. Note that in the 
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absence of flow, i// equals cod/c for every value of r. This figure characterizes the 
differences in results produced by the two methods. Theoretically, the results would be 
identical, but the different methods have different sources of error. The variation in 
results is particularly noticeable at small distances from the jet exit. 

The results presented herein have been calculated by using a relatively simple for- 
mulation for test purposes. The method is in no way restricted to the homogeneous con- 
vective wave equation or to the quasi -potential variable <£. Figure 8 illustrates the 



Figure 8.- Comparison of pressure formulation with qua si -potential formulation. 

results obtained with two alternate formulations. The label, M quasi-potential formulation 
with the extra term/’ refers to the equation which results when a correction term is added 
to the right-hand side of equation (1). Thus, 

L($)=i(v 2 u)|| (6) 

Pressure formulation implies use of the equation L(p) = 0 with the proper pressure 
information specified on the inner boundary. (See refs. 1, 2, and 6.) The comparison of 
pressure formulation with the two quasi- potential formulations were calculated at a flow 
Mach number of 0.025. Relative sound pressure level is plotted against 0. For this low 
flow case, the results are very much alike. 

Comparison of Numerical and Experimental Results 

Figure 9 is presented to illustrate the success of the current analysis for problems 
over a wide range of frequencies and flow Mach numbers. Numerical results from refer- 
ence 1 and experimental data of reference 4 are plotted as well as the numerical results 
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Present analysis 

Schubert (ref. 1) 

Experiment-Grande (ref. 4) 



(a) Effects of frequency. 



Figure 9*- Sound pressure level at r = lOOd and 0 = 0°. 


of the present method. Sound pressure level on the jet center line relative to SPL at 
6 = 90° is plotted against nondimensionalized frequency in figure 9(a) and against flow 
Mach number in figure 9(b). It is evident that the current analysis is in agreement with 
experimental findings. Moreover, the current method produces results which are even 
closer to experimental values than the results of reference 1. 

Providing more detailed information, the next two figures (figs. 10 and 11) show 
sample results at several flow velocities and acoustic source frequencies. In each graph, 
sound pressure level relative to SPL at 6 = 90° is plotted against 9 . These figures 
allow comparisons of the results of this paper with experimental data of reference 4 and 
the numerical results of reference 1. 


Figure 10 shows comparable results for the flow velocities Mach 0.3, Mach 0.5, and 
Mach 0.9. The radial distance from the source and the source frequency are lOOd and 
3000 Hz, respectively, in each of the three graphs. Notice that the correspondence between 
experimental data and the current numerical results is especially good at the lower Mach 
numbers. Also, notice the discrepancies between the numerical results of reference 1 
and experimental results at these same Mach numbers. 
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Figure 10.- Effects of jet exit Mach number on sound pressure level at r = lOOd. Source at 2d. 
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Numerical results of reference 5 are also presented in figure 10 for interest. Both 
the mathematical approach and the jet flow model used in reference 5 are significantly 
different from those of the current work. Thus, a fair comparison of one method with the 
other is not possible. 

The radial distance and flow velocity are constant (r = lOOd and Mj = 0.3) in fig- 
ure 11, while the frequency changes from 3000 to 5000 to 7000 Hz. In all cases, numerical 
results obtained by the present method represent an improvement over those obtained in 
reference 1 with respect to their agreement with experiment. 


Present analysis 

Schubert (ref. 1) 

Experiment -Grande (ref. 4) 



0 10 20 30 40 

e , deg 

(b) f = 5000 Hz. 

Figure 11.- Effects of frequency on sound pressure level at r = lOOd. Source at 2d. 
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Advantages of the Method 

It should be emphasized that the present method and the method of reference 1 use 
the same sound propagation model for testing purposes. Theoretically, they should pro- 
duce identical results. Differences arise when theory is translated into practical compu- 
ter programs. Iterative methods, such as the one described in reference 1, are often 
inaccurate because the iterative process must be terminated after a reasonable amount of 
computer time has elapsed. Furthermore, an iterative method is more susceptible to 
round-off errors which tend to multiply as computing time increases. Direct methods, 
such as the method presented herein, largely avoid these problems and generally produce 
more reliable results. 

The present method has other advantages besides accuracy. It is quite straight- 
forward and relatively easy to program. A typical program will execute in well under 
1 min of central memory time in a Control Data 6600 computer system. The computer 
storage space requirements are reasonable so that transfer of data to and from peripheral 
storage devices is unnecessary. 


CONCLUDING REMARKS 

Application of the method presented herein to the solution of sound propagation prob- 
lems has distinct advantages in that formulation of the problem in terms of a single slowly 
changing variable reduces the finite- difference approximation to a system of linear alge- 
braic equations. This system may then be solved by a standard library subroutine to give 
results which are accurate and yet economical to generate. Comparison with previous 
analysis and experimental data reveals that the present method yields more reliable results 
with decreased computer requirements. Finally, since the method is versatile and the 
programing task is relatively uncomplicated, the technique is a promising tool for future 
sound propagation studies. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
June 26, 1975 
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APPENDIX 


FORMULATION OF NUMERICAL SOLUTION 

The following steps were used to transform equation (1) for solution by numerical 
methods. First, equation (1) can be written in terms of nondimensional variables 

+ 2U*3>xT + tl^^xx) " A4> = 0 (Al) 

where T = tuj/d, C = c/uj, X = x/d, U = u/uj, and A is the normalized Laplacian 
operator. 

A transformation to spherical polar coordinates (R ,9,<p) causes the following changes 
in equation (Al): 


9 2 

3X 9T 


= cos 0 


9 2 

3T 9R 


- R 


-1 


sin 9 


9 2 

9T 90 


_a2_ 

3X 2 


= COS 


2 0 JL_ _ R-l sin 20 


9‘ 


-2 a 3 * 


9R^ 


3R 90 


+ R“* sin* 0 


90" 


+ R" 1 sin 2 0 + R 2 sin 20 ■— 

Furthermore, since it has been assumed that the jet is axisymmetric, ®<p(p = 0 an d the 
Laplacian operator A in spherical polar coordinates is defined by 

A — i 2 9 i 1 — i c °t 0 — 

9R.2 R 9R r2 802 r2 80 

Now, using equation (4) to specify <$, equation (Al) may be rewritten as 

a( K )z 


2 K 

I I <M R ’ 0 ) 

K= 0 1=0 


3R^ L ^ae^ 


= o 


(A2) 


where 


“0,0 = ■ M ^0 cos e) 2 + ^(l - 3 COS 2 0 ) - <// 0 2 

j[2MMjfi cos 0 + M 2 (l - 3 cos2 0)^ Q J 


■A 1 

R 2 
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q <1 o = - -^-(l - 3 cos2 o') + i^jQ MMjfi cos 9 + (l - M 2 cos 2 0 ^qJ 


cot 6 

ai > 1 


(l - 4M 2 sin 2 e) - i ^ 2M S1 ” - Mxp 0 cos 0) 


“ 2.0 cos2d ) 


a 0 = -Tr M 2 sin 20 
dfl R 2 


“ 2,2 = ^ ( X - M2 Sin20 ) 


The transformation into the (£, r\) coordinate system is defined by equation (5). 
applied to equation (A2) this transformation yields 


2 K 


y y « __aMz__ _ 

« s > iki 


where 


^0,0 = “o,0 

01.0 = R_1 (“l,0 - r_1 “2,o) 
^1,1 = B l(“l,l + B 2“2, 2^ 
^ 2,0 = R _ 2 “ 2,0 

02.1 = R ’ 1 ] B 1 “ 2,1 
02,2 = B l 2 “ 2,2 


and where 


Bl = 


q 2 sec 2 0 
(A + tan 2 0 


2 tan 0 (g^ - 1) 
cA + tan 2 0 


When 


(A3) 
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Finally, the set of difference equations on an (n - 2) X (m - 2) grid of the interior region 
of interest may be written 


I 


YK,L Z i+L-2,j-K+2 


= 0 


f (i = 2, 3, . . n - 1) 
\ (j = 2, 3, . . ., m - 1) 


(A4) 


where |y Kj can be defined as the 3x3 matrix 


02,1 

01,1 

02,1 

4HK 

2K k2 

4HK 

01,2 , 02,0 
2H h2 

P °’° K 2 H2 

01,0 02,0 
2H h 2 

02,1 

01,1 + 02,2 

,t hi 

4HK 

2K K 2 

4HK 


and where H = A£ and K= Arj. 

The additional equations needed to solve this system uniquely are the boundary con- 
ditions which follow. The inner boundary conditions are simply stated 


Z lfj = (AR)ij (j = 1, 2, . . m) 

Likewise, the discrete outer boundary conditions (refs. 1 and 2) are 

z n,j = " (!r‘i) Zn - 2 >i (j - 2, 3, . . m - 1) 

Lastly, the three -point symmetry conditions on the jet center line can be written 


Zi,l=|Zi,2-3Zi,3 


A 




(i = 2, 3, . . n) 


J i,m 


= ? Z 


i,m-l 3 z i,m-2 


(A5) 


(A 6 ) 


(A7) 
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